#### Least cost path calculation using passive telemetry datasets
## Written by VU 2018-02-12
## Last updated: 2020-03-20 (VU)
##
## Function to calculate least cost paths between detections from an ATTdata object (from VTrack)
## around landmasses. Uses 'osmdata' to construct cost raster if none is provided
## ************************************************************************************************************************************* ##
lcDistance <- function(ATTdata, cost = NULL, utm_epsg, ll_epsg = 4326, timestep = 60, h = 100, UDextent = 1, UDgrid = 100, cost.res = 100, directions = 16, ...){
####################################################################################
## ATTdata ATTdata object that consists of Station.Information, Tag.Detections and
## Tag.Metadata (output from setupData function in VTrack)
## cost cost raster, with land pixel values of 1000 and sea as 1, if none provided
## this is extracted and calculated using polygon downloaded from 'osmdata' package
## utm_epsg EPSG code for CRS object with projected coordinate reference system (units m)
## ll_epsg EPSG code for CRS object with unprojected coordinate reference system (units deg) [default 4326]
## timestep timestep for Center of activity positions (see COA() function in VTrack)
## h smoothing parameter for kernel density estimator in meters (default 100m)
## UDextent a value controlling the extent of the grid used for the estimation of kernel density (default = 4)
## UDgrid a number giving the size of the grid on which the kernel density should be estimated (in meters; default 100)
## cost.res resolution of cost layer (in degrees) calculated if not provided [default 100]
## directions number of directional axes to consider when calculating
## least cost path (options: 4, 8, 16 [default])
####################################################################################
## load required libraries and set up CRSs
sapply(c("lubridate","sf","dplyr","raster","gdistance","spatstat","maptools","rgeos","VTrack", "data.table", "lwgeom"), require, character.only=TRUE, warn.conflicts = F, quietly = T)
ll <- CRS(paste0("+init=epsg:", ll_epsg))
utm <- CRS(paste0("+init=epsg:", utm_epsg))
combdata <-
ATTdata$Tag.Detections %>%
left_join(ATTdata$Station.Information %>%
dplyr::select(Receiver, Station.Name, Station.Latitude, Station.Longitude),
by=c("Station.Name", "Receiver"))
statinfo <-
ATTdata$Station.Information %>%
st_as_sf(coords=c("Station.Longitude", "Station.Latitude"), crs = ll_epsg)
## Setup input tagdata; convert to spatial object and transform to projection
spdata_ll <-
combdata %>%
filter(!is.na(Station.Longitude)) %>%
mutate(step = 0:(nrow(.)-1)) %>%
st_as_sf(coords = c("Longitude", "Latitude"), crs = ll_epsg)
spdata_utm <-
spdata_ll %>%
st_transform(crs = utm_epsg)
## Extract landmass to calculate cost raster if not provided
if(is.null(cost)){
message("\n- No cost layer provided. Downloading coastline data from Open Street Map server")
require("osmdata")
tryCatch({
bb <- extent(statinfo) + 0.1
polydat <-
opq(bbox = bb[c(1,3,2,4)]) %>%
add_osm_feature(key = 'natural', value = 'coastline') %>%
osmdata_sf
islands <-
polydat$osm_polygons %>%
rowid_to_column() %>%
mutate(area = st_area(geometry)) %>%
dplyr::select(rowid, area)
coastline <- polydat$osm_lines %>% st_union %>% st_line_merge
pol <- st_as_sfc(st_bbox(coastline))
coastpoly <- lwgeom::st_split(st_geometry(pol), st_geometry(coastline))
pol_list <- list()
for(n in 1:length(coastpoly[[1]])){
pol_list[[n]] <- st_cast(coastpoly[[1]][n][[1]], 'POLYGON') %>% st_geometry() %>% st_sf()
st_crs(pol_list[[n]]) <- 4326
}
land <-
st_as_sf(data.table::rbindlist(pol_list)) %>%
rowid_to_column() %>%
mutate(area = st_area(geometry)) %>%
filter(area %in% max(area))
studypol <- st_as_sf(data.table::rbindlist(list(land, islands), use.names = T))
poly <-
studypol %>%
st_transform(utm_epsg) %>%
st_crop(extent(statinfo %>% st_transform(utm_epsg)) + 5000)
cost.ras<-rasterize(poly, raster(extent(statinfo %>% st_transform(utm_epsg)) + 5000, res = cost.res), 1000)
cost.ras[is.na(values(cost.ras))] <- 1
projection(cost.ras) <- utm
}, error=function(e){message("\nError: no land in sight!\nConsider adding your own cost layer")})
}else{
cost.in_utm <- raster::projectRaster(cost, crs=utm, method = "ngb")
cost.ras <- resample(cost.in_utm, raster(extent(cost.in_utm), res = cost.res), method = "ngb")
projection(cost.ras) <- utm
}
## Produce transition matrices, and correct for distortion
tryCatch({
message("\n- Calculating transition matrices for least cost path estimation")
trCost <- transition(1/cost.ras, mean, directions = directions)
trCost <- geoCorrection(trCost, type = "c")
}, error=function(e){message("\nError in calculating Transition layer")})
## Construct shortest path trajectories between sequence of detection steps
outdat <-
spdata_ll %>%
as_Spatial %>%
as_tibble %>%
rename(Longitude = coords.x1, Latitude = coords.x2) %>%
mutate(Distance_m = NA) %>%
dplyr::select(-step)
tryCatch({
message("- Constructing least cost path trajectories between consecutive detections")
traj<-list()
for(i in 1:max(spdata_utm$step)){
if(i %in% 1){pb <- txtProgressBar(min=1, max=max(spdata_utm$step), style=3)}
stepdat <- spdata_utm %>% filter(step %in% c(i-1, i))
origin <- stepdat %>% slice(1) %>% as_Spatial
goal <- stepdat %>% slice(n()) %>% as_Spatial
if (nrow(distinct(stepdat)) > 1){
traj[[i]]<-shortestPath(trCost, origin, goal, output = "SpatialLines")
outdat$Distance_m[i+1]<-as.numeric(costDistance(trCost, origin, goal))
} else {
traj[[i]]<- NULL
outdat$Distance_m[i+1]<-0
}
setTxtProgressBar(pb, i)
}
traj_clean <- Filter(Negate(is.null), traj)
trajectory_utm <- do.call(rbind, traj_clean)
trajectory_ll <- st_as_sf(trajectory_utm) %>% st_transform(ll_epsg)
}, error=function(e){message("\nError in calculating least cost path trajectories")})
## Calculate Kernel Density from trajectory
UDwin <- owin(xrange = (extent(spdata_utm) + (diff(extent(spdata_utm)[1:2])*UDextent))[1:2],
yrange = (extent(spdata_utm) + (diff(extent(spdata_utm)[3:4])*UDextent))[3:4])
dimyx <- c(diff(UDwin$yrange)/UDgrid,
diff(UDwin$xrange)/UDgrid)
# point UD
tryCatch({
message("\n- Estimating Kernel Density from COA positions")
coa.data <-
COA(ATTdata, timestep = timestep) %>%
st_as_sf(coords=c("Longitude.coa","Latitude.coa"), crs=4326, remove = F) %>%
st_transform(utm_epsg) %>% as_Spatial() %>% as_tibble() %>%
rename(X = coords.x1, Y = coords.x2)
suppressWarnings(pt.ppp <- ppp(x = coa.data$X, y = coa.data$Y, window = UDwin))
suppressWarnings(UDras.pt <- raster(density(pt.ppp, sigma = h, method = "C", dimyx = dimyx)))
values(UDras.pt) <- abs(values(UDras.pt)/max(values(UDras.pt)) - 1) * 100
projection(UDras.pt) <- utm
}, error=function(e){message("Error in estimating UD from COA positions")})
# trajectory UD
tryCatch({
message("- Estimating Kernel Density from least cost path trajectory")
traj.psp <- as.psp(trajectory_utm, window = UDwin)
UDras.traj <- raster(density(traj.psp, sigma = h, method = "C", dimyx = dimyx))
values(UDras.traj) <- abs(values(UDras.traj)/max(values(UDras.traj)) - 1) * 100
projection(UDras.traj) <- utm
}, error=function(e){message("Error in estimating UD for trajectories")})
# Merge trajectory UD with point UD and remove land areas
tryCatch({
UDras <- mean(UDras.traj, UDras.pt)
mask.ras <- resample(cost.ras, UDras, method = "ngb")
values(mask.ras)[values(mask.ras) > 1] <- NA
UDras.m_utm <- raster::mask(UDras, mask = mask.ras)
UDras.m_ll <- projectRaster(UDras.m_utm, crs = ll)
}, error=function(e){message("Error in estimating UD")})
# calculate UD areas
tryCatch({
message("- Estimating UD area for 50% and 95% contours")
UD50 <- UDras.m_utm < 50; UD50[values(UD50) %in% 0] <- NA
UD50.area <- gArea(rasterToPolygons(UD50, dissolve=T))
UD95 <- UDras.m_utm < 95; UD95[values(UD95) %in% 0] <- NA
UD95.area <- gArea(rasterToPolygons(UD95, dissolve=T))
}, error=function(e){message("Error in estimating UD areas")})
## return list output with step distances and spatial trajectory file
out<-list(tagdata = outdat,
coa.data = coa.data,
kernel.areas = data.frame(UD50_m2 = UD50.area, UD95_m2 = UD95.area),
lc.traj = trajectory_ll,
UD.raster = UDras.m_ll,
cost.raster = projectRaster(cost.ras, crs = ll, method = "ngb"))
return(out)
}
## ************************************************************************************************************************************* ##
## ************************************************************************************************************************************* ##
### Function to plot output on leaflet
plot.lcUD <- function(UDobj, ...){
library(mapview, warn.conflicts = F, quietly = T)
library(sf, warn.conflicts = F, quietly = T)
library(raster, warn.conflicts = F, quietly = T)
spdat <-
UDobj$tagdata %>%
st_as_sf(coords = c("Longitude", "Latitude"), crs = 4326) %>%
group_by(Station.Name) %>%
summarise(`Number of Detections` = n())
UDras <- UDobj$UD.raster
UDras[raster::values(UDras) > 99] <- NA
# coadat <-
# UDobj$coa.data %>%
# st_as_sf(coords = c("Longitude.coa", "Latitude.coa"), crs = 4326)
m <-
mapview(UDras, layer = "Kernel Utilisation Distribution", homebutton = F, na.color = "transparent", ...) +
# mapview(coadat, layer = "Center of Activity positions", col.regions = "white", alpha = 0, cex = 2, legend = F, homebutton = F) +
mapview(spdat, alpha = 0, layer = "Detection data", col.regions = "red", cex = "Number of Detections", legend = F, homebutton = F) +
mapview(st_geometry(UDobj$lc.traj), layer = "Least Cost Path Trajectories", legend = F, homebutton = F)
return(m)
}
## ************************************************************************************************************************************* ##
## ************************************************************************************************************************************* ##
## Function to generate least cost trajectories
lcTraj <- function(ATTdata, cost = NULL, utm_epsg, ll_epsg = 4326, cost.res = 50, directions = 16, ...){
## load required libraries and set up CRSs
sapply(c("lubridate","sf","dplyr","raster","gdistance","spatstat","maptools","rgeos","VTrack", "data.table", "lwgeom"), require,
character.only=TRUE, warn.conflicts = F, quietly = T)
ll <- CRS(paste0("+init=epsg:", ll_epsg))
utm <- CRS(paste0("+init=epsg:", utm_epsg))
combdata <-
ATTdata$Tag.Detections %>%
dplyr::select(-c(Longitude, Latitude)) %>%
left_join(ATTdata$Station.Information %>%
dplyr::select(Receiver, Station.Name, Station.Latitude, Station.Longitude),
by=c("Station.Name", "Receiver")) %>%
rename(Longitude = Station.Longitude,
Latitude = Station.Latitude)
statinfo <-
ATTdata$Station.Information %>%
st_as_sf(coords=c("Station.Longitude", "Station.Latitude"), crs = ll_epsg)
## Setup input tagdata; convert to spatial object and transform to projection
spdata_ll <-
combdata %>%
mutate(step = 0:(nrow(.)-1)) %>%
st_as_sf(coords = c("Longitude", "Latitude"), crs = ll_epsg)
spdata_utm <-
spdata_ll %>%
st_transform(crs = utm_epsg)
## Extract landmass to calculate cost raster if not provided
if(is.null(cost)){
message("- No cost layer provided. Downloading coastline data from Open Street Map server")
require("osmdata")
tryCatch({
bb <- extent(statinfo) + 0.1
polydat <-
opq(bbox = bb[c(1,3,2,4)]) %>%
add_osm_feature(key = 'natural', value = 'coastline') %>%
osmdata_sf
islands <-
polydat$osm_polygons %>%
rowid_to_column() %>%
mutate(area = st_area(geometry)) %>%
dplyr::select(rowid, area)
coastline <- polydat$osm_lines %>% st_union %>% st_line_merge
pol <- st_as_sfc(st_bbox(coastline))
coastpoly <- lwgeom::st_split(st_geometry(pol), st_geometry(coastline))
pol_list <- list()
for(n in 1:length(coastpoly[[1]])){
pol_list[[n]] <- st_cast(coastpoly[[1]][n][[1]], 'POLYGON') %>% st_geometry() %>% st_sf()
st_crs(pol_list[[n]]) <- 4326
}
land <-
st_as_sf(data.table::rbindlist(pol_list)) %>%
rowid_to_column() %>%
mutate(area = st_area(geometry)) %>%
filter(area %in% max(area))
studypol <- st_as_sf(data.table::rbindlist(list(land, islands), use.names = T))
poly <-
studypol %>%
st_transform(utm_epsg) %>%
st_crop(extent(statinfo %>% st_transform(utm_epsg)) + 5000)
cost.ras<-rasterize(poly, raster(extent(statinfo %>% st_transform(utm_epsg)) + 5000, res = cost.res), 1000)
cost.ras[is.na(values(cost.ras))] <- 1
projection(cost.ras) <- utm
}, error=function(e){message("Error: no land in sight!\nConsider adding your own cost layer")})
}else{
cost.in_utm <- raster::projectRaster(cost, crs = utm, method = "ngb")
cost.ras <- resample(cost.in_utm, raster(extent(cost.in_utm), res = cost.res), method = "ngb")
projection(cost.ras) <- utm
}
## Produce transition matrices, and correct for distortion
tryCatch({
message("- Calculating transition matrices for least cost path estimation")
trCost <- transition(1/cost.ras, mean, directions = directions)
trCost <- geoCorrection(trCost, type = "c")
}, error=function(e){message("Error in calculating Transition layer")})
## Construct shortest path trajectories between sequence of detection steps
outdat <-
spdata_ll %>%
as_Spatial %>%
as_tibble %>%
rename(Longitude = coords.x1, Latitude = coords.x2) %>%
mutate(Distance_m = NA) %>%
dplyr::select(-step)
tryCatch({
message("- Constructing least cost path trajectories between consecutive detections")
traj<-list()
for(i in 1:max(spdata_utm$step)){
if(i %in% 1){pb <- txtProgressBar(min=1, max=max(spdata_utm$step), style=3)}
stepdat <- spdata_utm %>% filter(step %in% c(i-1, i))
origin <- stepdat %>% slice(1) %>% as_Spatial
goal <- stepdat %>% slice(n()) %>% as_Spatial
if (nrow(distinct(stepdat)) > 1){
traj[[i]]<-shortestPath(trCost, origin, goal, output = "SpatialLines")
outdat$Distance_m[i+1]<-as.numeric(costDistance(trCost, origin, goal))
} else {
traj[[i]]<- NULL
outdat$Distance_m[i+1]<-0
}
setTxtProgressBar(pb, i)
}
traj_clean <- Filter(Negate(is.null), traj)
trajectory_utm <- do.call(rbind, traj_clean)
trajectory_ll <- st_as_sf(trajectory_utm) %>% st_transform(ll_epsg)
}, error=function(e){message("\nError in calculating least cost path trajectories")})
## return list output with step distances and spatial trajectory file
out<-list(tagdata = outdat,
lc.traj = trajectory_ll,
cost.raster = projectRaster(cost.ras, crs = ll, method = "ngb"))
return(out)
}
## ************************************************************************************************************************************* ##
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.